%     VIRTUAL UD - программа для получения эффективных свойств однонаправленного жгута
%     Copyright (C) 2020  Михаил Козлов
%     Сайт проекта: <http://www.fastenerdoctor.ru/>
%     Адрес электронной почты для связи с разработчиком:
%     <doctor.fastener@gmail.com>
% 
%     Эта программа является свободным программным обеспечением: Вы можете
%     распространять ее и/или изменять, соблюдая условия Универсальной общественной
%     лицензии GNU, опубликованной Фондом свободного программного обеспечения,
%     либо редакции 3, либо (на Ваше усмотрение) любой редакции, выпущенной позже.
% 
%     Эта программа распространяется в надежде на то, что она окажется полезной,
%     но БЕЗ КАКИХ-ЛИБО ГАРАНТИЙ, включая даже подразумеваемую гарантию ТОВАРНОЙ ПРИГОДНОСТИ,
%     либо ПРИГОДНОСТИ ДЛЯ ОПРЕДЕЛЕННЫХ ЦЕЛЕЙ. Ознакомьтесь с Универсальной общественной
%     лицензией GNU для получения более подробной информации.
% 
%     Вы должны были получить копию Универсальной общественной лицензии GNU
%     вместе с этой программой. Если Вы ее не получили, то перейдите по адресу:
%     <http://www.gnu.org/licenses/>.
% 
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
% 
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
%     GNU General Public License for more details.
% 
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <https://www.gnu.org/licenses/>.

function [NODES, ELEMENTS, RVE_Coordinates, Node_Merge_Tolerance] = Generate_Mesh(Fiber_Packing, Points, RVE_Thickness_Model, RVE_Thickness_Value, Element_Order, Mesh_Seed_Parameter, N_Elements_Thickness, Element_Critical_Ratio, RVE_Array_1, RVE_Array_2, RVE_Array_3, Analysis_Type)

global ErrCode;

global Merge_Tolerance_Multiplier;
global Auto_Thickness_Multiplier;

NODES = [];
ELEMENTS = [];
RVE_Coordinates = [];
Node_Merge_Tolerance = [];

fprintf (1, 'Построение сетки...                                     ');

switch Fiber_Packing
    
    case 'HEXAGONAL'
        
        [NODES_2D, ELEMENTS_2D, Min_Element_Edge, Max_Element_Edge, N_Layers_Suggestion] = Mesher_2D_Hexagonal(Points, Mesh_Seed_Parameter, Element_Critical_Ratio);

        if (ErrCode < 0)
            return;
        end
        
    case 'SQUARE'
        
        [NODES_2D, ELEMENTS_2D, Min_Element_Edge, Max_Element_Edge, N_Layers_Suggestion] = Mesher_2D_Square(Points, Mesh_Seed_Parameter, Element_Critical_Ratio);

        if (ErrCode < 0)
            return;
        end
        
    otherwise
        
        fprintf (1, '\nОшибка функции Generate_Mesh: недопустимое значение параметра Fiber_Packing !\n');
        ErrCode = -1;
        return;
    
end

% Автоматическое вычисление толщины ячейки
if (strcmp(RVE_Thickness_Model, 'AUTO'))
    
    switch Fiber_Packing
        
        case 'HEXAGONAL'
            
            RVE_Thickness_Value = norm(Points(18,:) - Points(55,:)) * Auto_Thickness_Multiplier;
            
        case 'SQUARE'
            
            RVE_Thickness_Value = norm(Points(10,:) - Points(29,:)) * Auto_Thickness_Multiplier;
        
    end
    
end

% Число элементов по высоте ячейки в случае автоматически вычисленной толщины
if (strcmp(RVE_Thickness_Model, 'AUTO'))
    N_Elements_Thickness = N_Layers_Suggestion;
end

% Уменьшение числа элементов по толщине ячейки в случае превышения максимально допустимого соотношения сторон элемента
Element_Ratio = Max_Element_Edge / (RVE_Thickness_Value/N_Elements_Thickness);

if (Element_Ratio > Element_Critical_Ratio)
    N_Elements_Thickness = floor(N_Elements_Thickness * Element_Critical_Ratio / Element_Ratio);
end

if (N_Elements_Thickness == 0)
    N_Elements_Thickness = 1;
end

% Увеличение числа элементов по толщине ячейки в случае превышения максимально допустимого соотношения сторон элемента
Element_Ratio = RVE_Thickness_Value/N_Elements_Thickness/Min_Element_Edge;

if (Element_Ratio > Element_Critical_Ratio)
    N_Elements_Thickness = ceil(N_Elements_Thickness * Element_Ratio / Element_Critical_Ratio);
end

% В механическом расчете должно быть не менее двух элементов по высоте, иначе не будет узлов для некоторых сэтов
if (strcmp(Analysis_Type, 'MECHANICAL') && (N_Elements_Thickness == 1))
    N_Elements_Thickness = 2;
end

% Обновление минимальной длины ребра элемента
Min_Element_Edge = min(Min_Element_Edge, RVE_Thickness_Value / N_Elements_Thickness);

% Допуск для слияния узлов
Node_Merge_Tolerance = Min_Element_Edge * Merge_Tolerance_Multiplier;

% Число узлов плоской сетки
N_Nodes_2D = size(NODES_2D,1);

% Выделение памяти под узлы объемной сетки
NODES = zeros((N_Elements_Thickness + 1) * N_Nodes_2D, 4);

% Перестановка координат в массиве узлов плоской сетки
NODES_2D(:,4) = NODES_2D(:,3);
NODES_2D(:,3) = NODES_2D(:,2);
NODES_2D(:,2) = zeros(N_Nodes_2D,1);

% Заполнение массива узлов объемной сетки
for i=0:N_Elements_Thickness
    
    NODES(i*N_Nodes_2D + 1:(i+1)*N_Nodes_2D,:) = NODES_2D(:,:);
    
    if (i > 0)
        
        NODES(i*N_Nodes_2D + 1:(i+1)*N_Nodes_2D, 1) = NODES((i-1)*N_Nodes_2D + 1:i*N_Nodes_2D, 1) + ones(N_Nodes_2D,1) * N_Nodes_2D;
        NODES(i*N_Nodes_2D + 1:(i+1)*N_Nodes_2D, 2) = ones(N_Nodes_2D,1) * i/N_Elements_Thickness * RVE_Thickness_Value;
        
    end
    
end

% Число элементов плоской сетки
N_Elements_2D = size(ELEMENTS_2D,1);

% Выделение памяти под элементы объемной сетки
ELEMENTS = zeros(N_Elements_Thickness * N_Elements_2D, 12);

% Заполнение первого слоя объемной сетки
for i=1:N_Elements_2D
            
    switch ELEMENTS_2D(i,3)

        case 3

            ELEMENTS(i,1) = ELEMENTS_2D(i,1);
            ELEMENTS(i,2) = 3;
            ELEMENTS(i,3) = 6;
            ELEMENTS(i,4) = ELEMENTS_2D(i,4);

            Element_Normal = cross(NODES_2D(ELEMENTS_2D(i,6),2:4) - NODES_2D(ELEMENTS_2D(i,5),2:4), NODES_2D(ELEMENTS_2D(i,7),2:4) - NODES_2D(ELEMENTS_2D(i,5),2:4));

            if (dot(Element_Normal, [1,0,0]) > 0)

                ELEMENTS(i,5) = ELEMENTS_2D(i,5);
                ELEMENTS(i,6) = ELEMENTS_2D(i,6);
                ELEMENTS(i,7) = ELEMENTS_2D(i,7);
                ELEMENTS(i,8) = ELEMENTS_2D(i,5) + N_Nodes_2D;
                ELEMENTS(i,9) = ELEMENTS_2D(i,6) + N_Nodes_2D;
                ELEMENTS(i,10) = ELEMENTS_2D(i,7) + N_Nodes_2D;

            else

                ELEMENTS(i,5) = ELEMENTS_2D(i,7);
                ELEMENTS(i,6) = ELEMENTS_2D(i,6);
                ELEMENTS(i,7) = ELEMENTS_2D(i,5);
                ELEMENTS(i,8) = ELEMENTS_2D(i,7) + N_Nodes_2D;
                ELEMENTS(i,9) = ELEMENTS_2D(i,6) + N_Nodes_2D;
                ELEMENTS(i,10) = ELEMENTS_2D(i,5) + N_Nodes_2D;

            end

        case 4

            ELEMENTS(i,1) = ELEMENTS_2D(i,1);
            ELEMENTS(i,2) = 3;
            ELEMENTS(i,3) = 8;
            ELEMENTS(i,4) = ELEMENTS_2D(i,4);

            Element_Normal = cross(NODES_2D(ELEMENTS_2D(i,6),2:4) - NODES_2D(ELEMENTS_2D(i,5),2:4), NODES_2D(ELEMENTS_2D(i,8),2:4) - NODES_2D(ELEMENTS_2D(i,5),2:4));

            if (dot(Element_Normal, [1,0,0]) > 0)

                ELEMENTS(i,5) = ELEMENTS_2D(i,5);
                ELEMENTS(i,6) = ELEMENTS_2D(i,6);
                ELEMENTS(i,7) = ELEMENTS_2D(i,7);
                ELEMENTS(i,8) = ELEMENTS_2D(i,8);
                ELEMENTS(i,9) = ELEMENTS_2D(i,5) + N_Nodes_2D;
                ELEMENTS(i,10) = ELEMENTS_2D(i,6) + N_Nodes_2D;
                ELEMENTS(i,11) = ELEMENTS_2D(i,7) + N_Nodes_2D;
                ELEMENTS(i,12) = ELEMENTS_2D(i,8) + N_Nodes_2D;

            else

                ELEMENTS(i,5) = ELEMENTS_2D(i,8);
                ELEMENTS(i,6) = ELEMENTS_2D(i,7);
                ELEMENTS(i,7) = ELEMENTS_2D(i,6);
                ELEMENTS(i,8) = ELEMENTS_2D(i,5);
                ELEMENTS(i,9) = ELEMENTS_2D(i,8) + N_Nodes_2D;
                ELEMENTS(i,10) = ELEMENTS_2D(i,7) + N_Nodes_2D;
                ELEMENTS(i,11) = ELEMENTS_2D(i,6) + N_Nodes_2D;
                ELEMENTS(i,12) = ELEMENTS_2D(i,5) + N_Nodes_2D;

            end

    end
    
end

% Заполнение остальных слоев объемной сетки
Mask = ELEMENTS(1 : N_Elements_2D, 5:12) > 0;
Mask = double(Mask);
    
for i=1:N_Elements_Thickness-1
    
    ELEMENTS(i*N_Elements_2D + 1 : (i+1)*N_Elements_2D, 1) = ELEMENTS((i-1)*N_Elements_2D + 1 : i*N_Elements_2D, 1) + ones(N_Elements_2D,1)*N_Elements_2D;
    ELEMENTS(i*N_Elements_2D + 1 : (i+1)*N_Elements_2D, 2:4) = ELEMENTS((i-1)*N_Elements_2D + 1 : i*N_Elements_2D, 2:4);
    ELEMENTS(i*N_Elements_2D + 1 : (i+1)*N_Elements_2D, 5:12) = ELEMENTS((i-1)*N_Elements_2D + 1 : i*N_Elements_2D, 5:12) + Mask*N_Nodes_2D;
    
end

% Габариты ячейки
RVE_Coordinates = zeros(3,2);

switch Fiber_Packing

    case 'HEXAGONAL'

        RVE_Coordinates(1,1) = 0;
        RVE_Coordinates(1,2) = RVE_Thickness_Value;
        RVE_Coordinates(2,1) = Points(11,1);
        RVE_Coordinates(2,2) = Points(32,1);
        RVE_Coordinates(3,1) = Points(35,2);
        RVE_Coordinates(3,2) = Points(8,2);

    case 'SQUARE'

        RVE_Coordinates(1,1) = 0;
        RVE_Coordinates(1,2) = RVE_Thickness_Value;
        RVE_Coordinates(2,1) = Points(3,1);
        RVE_Coordinates(2,2) = Points(18,1);
        RVE_Coordinates(3,1) = Points(26,2);
        RVE_Coordinates(3,2) = Points(2,2);

end

fprintf (1, 'OK!\n');

% Конвертация элементов в элементы второго порядка
if (strcmp(Element_Order, 'QUADRATIC'))
    
    fprintf (1, 'Конвертация элементов в квадратичные...                 ');
    
    % Уменьшение допуска для слияния узлов
    Node_Merge_Tolerance = Node_Merge_Tolerance/2;
    
    [NODES, ELEMENTS] = Convert_Mesh_To_Quad(NODES, ELEMENTS, Node_Merge_Tolerance);
    
    % Проецирование срединных узлов на границу волокон
    switch Fiber_Packing
        
        case 'HEXAGONAL'
            
            R = Points(5,2);
            Shift_Y = -Points(11,1);
            Shift_Z = Points(8,2);
            
            for i=1:size(NODES,1)
                
                % Центральное волокно
            
                for j=1:3*Mesh_Seed_Parameter
                    
                    Midnode_Coord = ([cosd(30/Mesh_Seed_Parameter * j), sind(30/Mesh_Seed_Parameter * j)] + [cosd(30/Mesh_Seed_Parameter * (j-1)), sind(30/Mesh_Seed_Parameter * (j-1))]) / 2 * R;
                
                    Point_Coord = [cosd(30/Mesh_Seed_Parameter * j - 15/Mesh_Seed_Parameter), sind(30/Mesh_Seed_Parameter * j - 15/Mesh_Seed_Parameter)] * R;
                    
                    if ((abs(NODES(i,3) - Midnode_Coord(1)) < Node_Merge_Tolerance) && (abs(NODES(i,4) - Midnode_Coord(2)) < Node_Merge_Tolerance))
                        NODES(i,3) = Point_Coord(1);
                        NODES(i,4) = Point_Coord(2);
                    end
                
                end
                
                for j=1:3*Mesh_Seed_Parameter
                    
                    Midnode_Coord = ([-cosd(30/Mesh_Seed_Parameter * j), -sind(30/Mesh_Seed_Parameter * j)] + [-cosd(30/Mesh_Seed_Parameter * (j-1)), -sind(30/Mesh_Seed_Parameter * (j-1))]) / 2 * R;
                
                    Point_Coord = [-cosd(30/Mesh_Seed_Parameter * j - 15/Mesh_Seed_Parameter), -sind(30/Mesh_Seed_Parameter * j - 15/Mesh_Seed_Parameter)] * R;
                    
                    if ((abs(NODES(i,3) - Midnode_Coord(1)) < Node_Merge_Tolerance) && (abs(NODES(i,4) - Midnode_Coord(2)) < Node_Merge_Tolerance))
                        NODES(i,3) = Point_Coord(1);
                        NODES(i,4) = Point_Coord(2);
                    end
                
                end
                
                for j=1:3*Mesh_Seed_Parameter
                    
                    Midnode_Coord = ([-cosd(30/Mesh_Seed_Parameter * j), sind(30/Mesh_Seed_Parameter * j)] + [-cosd(30/Mesh_Seed_Parameter * (j-1)), sind(30/Mesh_Seed_Parameter * (j-1))]) / 2 * R;
                
                    Point_Coord = [-cosd(30/Mesh_Seed_Parameter * j - 15/Mesh_Seed_Parameter), sind(30/Mesh_Seed_Parameter * j - 15/Mesh_Seed_Parameter)] * R;
                    
                    if ((abs(NODES(i,3) - Midnode_Coord(1)) < Node_Merge_Tolerance) && (abs(NODES(i,4) - Midnode_Coord(2)) < Node_Merge_Tolerance))
                        NODES(i,3) = Point_Coord(1);
                        NODES(i,4) = Point_Coord(2);
                    end
                
                end
                
                for j=1:3*Mesh_Seed_Parameter
                    
                    Midnode_Coord = ([cosd(30/Mesh_Seed_Parameter * j), -sind(30/Mesh_Seed_Parameter * j)] + [cosd(30/Mesh_Seed_Parameter * (j-1)), -sind(30/Mesh_Seed_Parameter * (j-1))]) / 2 * R;
                
                    Point_Coord = [cosd(30/Mesh_Seed_Parameter * j - 15/Mesh_Seed_Parameter), -sind(30/Mesh_Seed_Parameter * j - 15/Mesh_Seed_Parameter)] * R;
                    
                    if ((abs(NODES(i,3) - Midnode_Coord(1)) < Node_Merge_Tolerance) && (abs(NODES(i,4) - Midnode_Coord(2)) < Node_Merge_Tolerance))
                        NODES(i,3) = Point_Coord(1);
                        NODES(i,4) = Point_Coord(2);
                    end
                
                end
                
                % Угловые волокна
                
                for j=1:3*Mesh_Seed_Parameter
                    
                    Midnode_Coord = ([cosd(30/Mesh_Seed_Parameter * j), sind(30/Mesh_Seed_Parameter * j)] + [cosd(30/Mesh_Seed_Parameter * (j-1)), sind(30/Mesh_Seed_Parameter * (j-1))]) / 2 * R + [-Shift_Y, -Shift_Z];
                
                    Point_Coord = [cosd(30/Mesh_Seed_Parameter * j - 15/Mesh_Seed_Parameter), sind(30/Mesh_Seed_Parameter * j - 15/Mesh_Seed_Parameter)] * R + [-Shift_Y, -Shift_Z];
                    
                    if ((abs(NODES(i,3) - Midnode_Coord(1)) < Node_Merge_Tolerance) && (abs(NODES(i,4) - Midnode_Coord(2)) < Node_Merge_Tolerance))
                        NODES(i,3) = Point_Coord(1);
                        NODES(i,4) = Point_Coord(2);
                    end
                
                end
                
                for j=1:3*Mesh_Seed_Parameter
                    
                    Midnode_Coord = ([-cosd(30/Mesh_Seed_Parameter * j), -sind(30/Mesh_Seed_Parameter * j)] + [-cosd(30/Mesh_Seed_Parameter * (j-1)), -sind(30/Mesh_Seed_Parameter * (j-1))]) / 2 * R + [Shift_Y, Shift_Z];
                
                    Point_Coord = [-cosd(30/Mesh_Seed_Parameter * j - 15/Mesh_Seed_Parameter), -sind(30/Mesh_Seed_Parameter * j - 15/Mesh_Seed_Parameter)] * R + [Shift_Y, Shift_Z];
                    
                    if ((abs(NODES(i,3) - Midnode_Coord(1)) < Node_Merge_Tolerance) && (abs(NODES(i,4) - Midnode_Coord(2)) < Node_Merge_Tolerance))
                        NODES(i,3) = Point_Coord(1);
                        NODES(i,4) = Point_Coord(2);
                    end
                
                end
                
                for j=1:3*Mesh_Seed_Parameter
                    
                    Midnode_Coord = ([-cosd(30/Mesh_Seed_Parameter * j), sind(30/Mesh_Seed_Parameter * j)] + [-cosd(30/Mesh_Seed_Parameter * (j-1)), sind(30/Mesh_Seed_Parameter * (j-1))]) / 2 * R + [Shift_Y, -Shift_Z];
                
                    Point_Coord = [-cosd(30/Mesh_Seed_Parameter * j - 15/Mesh_Seed_Parameter), sind(30/Mesh_Seed_Parameter * j - 15/Mesh_Seed_Parameter)] * R + [Shift_Y, -Shift_Z];
                    
                    if ((abs(NODES(i,3) - Midnode_Coord(1)) < Node_Merge_Tolerance) && (abs(NODES(i,4) - Midnode_Coord(2)) < Node_Merge_Tolerance))
                        NODES(i,3) = Point_Coord(1);
                        NODES(i,4) = Point_Coord(2);
                    end
                
                end
                
                for j=1:3*Mesh_Seed_Parameter
                    
                    Midnode_Coord = ([cosd(30/Mesh_Seed_Parameter * j), -sind(30/Mesh_Seed_Parameter * j)] + [cosd(30/Mesh_Seed_Parameter * (j-1)), -sind(30/Mesh_Seed_Parameter * (j-1))]) / 2 * R + [-Shift_Y, Shift_Z];
                
                    Point_Coord = [cosd(30/Mesh_Seed_Parameter * j - 15/Mesh_Seed_Parameter), -sind(30/Mesh_Seed_Parameter * j - 15/Mesh_Seed_Parameter)] * R + [-Shift_Y, Shift_Z];
                    
                    if ((abs(NODES(i,3) - Midnode_Coord(1)) < Node_Merge_Tolerance) && (abs(NODES(i,4) - Midnode_Coord(2)) < Node_Merge_Tolerance))
                        NODES(i,3) = Point_Coord(1);
                        NODES(i,4) = Point_Coord(2);
                    end
                
                end
            
            end
            
        case 'SQUARE'
            
            R = Points(2,2) + Points(4,1);
            Shift_Y = Points(2,2);
            Shift_Z = Points(2,2);
            
            for i=1:size(NODES,1)
                
                % Угловые волокна
                
                for j=1:2*Mesh_Seed_Parameter
                    
                    Midnode_Coord = ([cosd(45/Mesh_Seed_Parameter * j), sind(45/Mesh_Seed_Parameter * j)] + [cosd(45/Mesh_Seed_Parameter * (j-1)), sind(45/Mesh_Seed_Parameter * (j-1))]) / 2 * R + [-Shift_Y, -Shift_Z];
                
                    Point_Coord = [cosd(45/Mesh_Seed_Parameter * j - 22.5/Mesh_Seed_Parameter), sind(45/Mesh_Seed_Parameter * j - 22.5/Mesh_Seed_Parameter)] * R + [-Shift_Y, -Shift_Z];
                    
                    if ((abs(NODES(i,3) - Midnode_Coord(1)) < Node_Merge_Tolerance) && (abs(NODES(i,4) - Midnode_Coord(2)) < Node_Merge_Tolerance))
                        NODES(i,3) = Point_Coord(1);
                        NODES(i,4) = Point_Coord(2);
                    end
                
                end
                
                for j=1:2*Mesh_Seed_Parameter
                    
                    Midnode_Coord = ([-cosd(45/Mesh_Seed_Parameter * j), -sind(45/Mesh_Seed_Parameter * j)] + [-cosd(45/Mesh_Seed_Parameter * (j-1)), -sind(45/Mesh_Seed_Parameter * (j-1))]) / 2 * R + [Shift_Y, Shift_Z];
                
                    Point_Coord = [-cosd(45/Mesh_Seed_Parameter * j - 22.5/Mesh_Seed_Parameter), -sind(45/Mesh_Seed_Parameter * j - 22.5/Mesh_Seed_Parameter)] * R + [Shift_Y, Shift_Z];
                    
                    if ((abs(NODES(i,3) - Midnode_Coord(1)) < Node_Merge_Tolerance) && (abs(NODES(i,4) - Midnode_Coord(2)) < Node_Merge_Tolerance))
                        NODES(i,3) = Point_Coord(1);
                        NODES(i,4) = Point_Coord(2);
                    end
                
                end
                
                for j=1:2*Mesh_Seed_Parameter
                    
                    Midnode_Coord = ([-cosd(45/Mesh_Seed_Parameter * j), sind(45/Mesh_Seed_Parameter * j)] + [-cosd(45/Mesh_Seed_Parameter * (j-1)), sind(45/Mesh_Seed_Parameter * (j-1))]) / 2 * R + [Shift_Y, -Shift_Z];
                
                    Point_Coord = [-cosd(45/Mesh_Seed_Parameter * j - 22.5/Mesh_Seed_Parameter), sind(45/Mesh_Seed_Parameter * j - 22.5/Mesh_Seed_Parameter)] * R + [Shift_Y, -Shift_Z];
                    
                    if ((abs(NODES(i,3) - Midnode_Coord(1)) < Node_Merge_Tolerance) && (abs(NODES(i,4) - Midnode_Coord(2)) < Node_Merge_Tolerance))
                        NODES(i,3) = Point_Coord(1);
                        NODES(i,4) = Point_Coord(2);
                    end
                
                end
                
                for j=1:2*Mesh_Seed_Parameter
                    
                    Midnode_Coord = ([cosd(45/Mesh_Seed_Parameter * j), -sind(45/Mesh_Seed_Parameter * j)] + [cosd(45/Mesh_Seed_Parameter * (j-1)), -sind(45/Mesh_Seed_Parameter * (j-1))]) / 2 * R + [-Shift_Y, Shift_Z];
                
                    Point_Coord = [cosd(45/Mesh_Seed_Parameter * j - 22.5/Mesh_Seed_Parameter), -sind(45/Mesh_Seed_Parameter * j - 22.5/Mesh_Seed_Parameter)] * R + [-Shift_Y, Shift_Z];
                    
                    if ((abs(NODES(i,3) - Midnode_Coord(1)) < Node_Merge_Tolerance) && (abs(NODES(i,4) - Midnode_Coord(2)) < Node_Merge_Tolerance))
                        NODES(i,3) = Point_Coord(1);
                        NODES(i,4) = Point_Coord(2);
                    end
                
                end
            
            end
            
        otherwise
            
            fprintf (1, '\nОшибка функции Generate_Mesh: недопустимое значение параметра Fiber_Packing !\n');
            ErrCode = -1;
            return;
        
    end
    
    fprintf (1, 'OK!\n');
    
end

% Создание массива ячеек
if ((RVE_Array_1 > 1) || (RVE_Array_2 > 1) || (RVE_Array_3 > 1))
    
    fprintf (1, 'Создание массива ячеек...                               ');
    
    % Размеры одинарной ячейки периодичности
    RVE_Size = zeros(1,3);
    
    switch Fiber_Packing

        case 'HEXAGONAL'

            RVE_Size(1) = RVE_Thickness_Value;
            RVE_Size(2) = Points(32,1) - Points(11,1);
            RVE_Size(3) = Points(8,2) - Points(35,2);

        case 'SQUARE'

            RVE_Size(1) = RVE_Thickness_Value;
            RVE_Size(2) = Points(18,1) - Points(3,1);
            RVE_Size(3) = Points(2,2) - Points(26,2);

    end
    
    [NODES, ELEMENTS, RVE_Coordinates] = Create_RVE_Array(NODES, ELEMENTS, Node_Merge_Tolerance, RVE_Array_1, RVE_Array_2, RVE_Array_3, RVE_Size, RVE_Coordinates);
    
    fprintf (1, 'OK!\n');
    
end
